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Abstract 


This  paper  discusses  in  substantial  detail  the  design  principles 
and  structure  of  an  existing  Fortran  program  library  whose  primary  appli- 
cation is  to  solve  optimization  problems.  Such  a discussion  not  only  helps 
to  clarify  the  scope  of  application  for  potential  users  of  the  library,  but 
also  is  useful  for  workers  on  other  software  projects.  The  fundamental 
objectives  of  the  present  library  have  been  to  produce  sound,  careful  imple- 
mentations of  reliable  methods  that  represent  the  state  of  the  art  in  numeri- 
cal optimization.  The  general  implications  of  these  overall  design  aims 
are  presented,  as  well  as  specific  instances  of  the  results  of  decisions  to 
include  particular  desirable  features. 
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^ ^ The  evoIyinK  need  for  pro>;ram  l ibraries 

The  manner  in  which  computer  implementations  of  numerical  algorithms 
are  disseminated  has  changed  dramatically  since  the  early  days  of  auto- 
matic computation.  Initially ^ a user  who  required  a computer  routine  to 
solve  a particular  numerical  problem  would  typically  consult  research 
journals  that  might  contain  a theoretical  description  of  an  appropriate 
method,  and  then  write  his  own  code.  Unfortunately,  such  personalized 
implementations  were  subject  to  a significant  risk  of  unreliability, 
because  of  a lack  of  attention  to  details  of  programming  or  numerical 
analysis.  Furthermore,  as  the  complexity  of  numerical  methods  increased, 
it  became  impractical  for  individuals  to  write  their  own  versions  of  all 
necessary  algorithms,  even  given  the  will  to  do  so. 

Subsequently,  it  became  the  practice  among  some  authors  of  new  numeri- 
cal methods  to  publish  computer  programs  as  well  as  theoretical  descrip- 
tions. Although  this  development  was  a step  in  the  right  direction,  in 
many  ways  the  situation  was  even  more  complicated.  The  quality  of  the 
published  programs  varied  enormously,  and  there  was  little  uniformity  of 
standards  concerning  programming  structure  and  style.  In  addition,  the 
publi' hed  programs  were  often  written  only  to  test  the  proposed  method  on 
a small  set  of  problems  (usually  well-behaved),  and  displayed  serious  errors 
of  prograniming  as  well  as  the  inability  to  detect  and  recover  from  numerical 
difficulties.  Under  these  conditions,  the  user  was  obliged  to  undertake  a 
search  of  the  literature,  among  a large  collection  ol  published  routines, 
with  no  guidelines  to  assist  in  making  a good  choice. 
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The  inefficiency  inJ  confusion  resulting  from  an  uncontrolled  prolif- 


e.ration of  alternative  routines  led  to  an  awareness  of  tlie  need  for  program 
libraries . This  term  does  not  moan  simply  a collection  of  programs  from 
varied  sources,  written  in  isolation,  thai:  are  grouped  together  in  one  place 
with  some  uniform  documentation.  Rather,  a program  library  is  a sot  of  rou- 
tines that  are  conceived  and  written  within  a unified  framework,  to  be  avail- 
able to  a general  community  of  users.  Today  there  is  wide  recognition  of 
the  great  value  of  program  libraries  designed  to  solve  useful  classes  of 
numerical  problems,  since  the  latest  developments  in  mathematical  computing 
can  thereby  be  made  available  to  many  users. 

1 .2  Issues  in  the  design  of  a program  library 

Two  aspects  of  the  design  and  production  of  a program  library  will  be 
considered  here. 

First,  it  is  essential  that  each  individual  routine  contained  in  a 
program  library  can  be  classified  as  "good"  software.  This  requirement 
implies  the  need  for  analysis  of  the  purpose  of  every  routine,  since  the 
task  of  developing  a sound  and  careful  implementation  of  a numerical  method 
is  extremely  difficult  and  time-consuming,  even  for  experts  i see  the  Preface 
to  Wilkinson  and  Reinsch,  I9TI principles  upon  which  good  computer 
programs  for  numerical  methods  should  be  based  have  been  discussed  by 
numerous  authors  in  varying  contexts  (see,  tor  example,  Byrne  and  Hindmarsh, 
1975;  Cody,  1971,  ; Cox,  197!  ; do.  Boor,  1971a;  Dixon,  1 97)  ; Ford 

and  Hague,  197'*;  Hillstrom  et  al.,  107'  ; Rice,  1971;  Schonfelder,  197"; 
Shampine  and  Gordon,  I9''.  ; Smith  et  al.,  19*713,'.  Almost  without  exception, 
it  is  agreed  that  mathematical  software  should  be  designed  to  satisfy  certain 
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criteria  (see  Rice,  L *'(1^  lor  a widely  fjuoLed  list).  However,  it  is  further 
recognized  that  some  of  these  desirable  qualities  are  inherently  contra- 
dictory (see,  lor  example,  Cox,  hixon,  1 yY'. ; Hillstrom  et  al.,  l^'Y'  ; 

Kahaner,  l'V('l,  lO'flj;  l.yness  and  Kag.anove,  i'/ft;  Rice,  1971).  Consequently, 
the  creation  of  any  computer  program  must  include  decisions,  implicit  and 
explicit,  concerning  ttie  reiativ'c  weight  and  iiiiportance  to  be  assigned  to 
Che  possibly  conflicting  attributes. 

Second,  a uniform  global  design  should  be  imposed  on  the  entire  program 
library,  wiiose  definition  implies  the  existence  of  a colfective  purpose  that 
transcends  the  aims  of  any  particular  routine.  The  significant  assumption 
that  a program  library  sliouid  bo  available  and  useful  to  a general-user 
community  has  several  inimediate  consequences.  In  particular,  the  demands 
upon  the  library  routines  will  necessarily  reflect  diflering  levels  of 
interest  and  expertise,  and  hence  will  vary  widely  from  user  to  user.  For 
example,  a user  who  is  familiar  with  the  details  of  an  algorithm  may  wish 
the  corresponding  library  imj-ilementacion  to  contain  sufficient  flexibility 
to  be  "tuned"  to  a particular  frequently  occurring  problem;  a user  who 
simply  needs  Co  solve  a problem  one  time  only  may  want  to  be  able  to  use 
Che  library  routines  with  a minimum  of  time  and  effort,  without  being 
required  to  learn  any  details  of  tlie  underlying  algorithms  or  software. 
Similarly,  the  criteria  that  are  considered  desirable  in  the  library  rou- 
tines will  change  drastically  with  different  applications.  In  one  instance, 
a library  program  may  be  used  to  solve  a problem  for  which  the  cost  of  com- 
puter time  is  negligible  compared  to  the  implications  of  failing  to  solve 
the.  problem  or  of  finding  a poor  solution,  so  that  the  need  for  reliability 
dominates  all  other  criteria.  In  another  application  of  the  same  numerical 


problum,  however,  the  most  important  consideration  may  be  spcicd  of  execution, 
even  at  the  risk  of  inaccuracy  or  faiLurc!. 

Clearly,  a "perfect"  library  can  never  exist.  Any  particular  routine, 
or  the  entire  library,  will  inevitably  be  subject  to  criticism  from  some 
individuals  for  not  satisfying  their  ideal  requirements,  since  all  desirable 
features  cannot  be  adiieved  simultaneously.  Therefore,  the  design  decisions 
and  the  justification  of  those  decisions  should  be  specified  explicitly  for 
any  program  library.  A presentation  of  the  overall  philosophy  behind  a 
library  helps  to  clarify  the  scope  of  application  for  potential  users,  and 
allows  an  informed  judgment  of  how  the  available  routines  may  be  used  most 
effectively  for  any  given  problem.  In  addition,  such  a discussion  is  help- 
ful for  workers  on  other  software  projects  because  similar  difficulties  and 
questions  often  recur  in  efforts  to  create  program  libraries;  a discussion 
of  experience  and  alternatives  in  a particular  library  design  can  reveal 
non-obvious,  or  unfortunate,  consequences  of  certain  critical  decisions, 
which  might  not  otherwise  become  apparent  until  significant  effort  had  been 
expended , 

1 .5  Purpose  of  this  paper 

Descriptions  of  collections  of  computer  routines  have  been  published, 
in  varying  degrees  of  detail,  for  several  areas  of  numerical  analysis: 
quadrature  (see  de  Boor,  1971b);  data  fitting  (see  Cox,  l97li ) ; ordinary 
differential  equations  (see  Byrne  and  Hindmarsh,  I975;  Shampine  and  Gordon, 
^977);  nonlinear  partial  differential  equations  (see  Sincovec  and  Madsen, 
1975);  eigensystems  (sec  Smith  et  al.,  197^ib ) ; and  special  functions  (see 
Cody,  197'  )•  In  some  of  these  cases,  the  motivation  for  the  chosen  imple- 
mentation is  presented,  either  explicitly  or  by  Implication.  There  are 


also  examples  of  documentation  of  routines  for  optimization  (see,  for 
example,  Fletcher,  1972;  James  and  Roos , 19  ([3;  hand  and  Powell,  1973! 
Lasdon  et  al.,  1970;  Mylandcr,  Holmes,  and  McCormick,  197^^!  Rosen  and 
Wagner,  1975) • However,  none  of  these  latter  instances  includes  a dis- 
cussion of  the  reasons  for  the  given  design,  nor  qualifies  as  a "program 
library” . 

The  purpose  of  this  paper  is  to  discuss  in  substantial  detail  the 
design  principles  and  structure  of  an  existing  Fortran  library  whose 
primary  application  is  to  solve  optimization  problems.  The  detail  is 
included  to  avoid  general  statements  of  principles  with  which  everyone 
would  agree.  There  is  an  immense  gap  between  a theoretical  description 
of  intent  in  software  vrriting  — for  example,  "we  emphasize  modularity, 
good  structure,  and  efficiency"  — and  an  analysis  of  how  well  the  chosen 
design  in  fact  satisfies  the  proposed  criteria.  Only  the  latter  process 
determines  the  quality  3 resulting  software. 

This  paper  is  n o serve  as  documentation  for  the  library, 

but  rather  to  provi  rspective  on  the  complete  design;  however, 

specific  examples  will  always  be  cited  to  illustrate  relevant  points. 
Detailed  documentation  of  the  routines  is  published  separately  (see,  for 
example.  Gill  et  al.,  1975^^  3nd  the  appendix). 
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2 . Overall  Philosophy 

2 . 1 Statement  of  problem  and  objectives  of  library  design 

The  optimization  problems  to  be  considered  can  be  stated  as  follows: 

PI:  minimize  F(x),  x £ 

subject  to  c^(x)  - 0 , i --  l,2,...,m', 
c^(x)  ^ 0 ^ i = m''+l^...jm  , 

where  F and  given  real-valued  functions.  The  function  F(x) 

is  normally  termed  the  "objective  function,"  and  the  set  {c^}  is  the  set 
of  "constraint  functions." 

Certain  kinds  of  numerical  problems  have  the  property  that  the  algo- 
rithms to  solve  them  are  finite  decision  procedures,  or  "reliable  exact 
arithmetic  algorithms"  (see  Lyness  and  Kaganove,  19?6,  for  a further  dis- 
cussion of  ways  to  characterize  numerical  problems).  To  create  a good 
computer  implementation  of  such  an  algorithm  is  in  itself  a difficult  task, 
but  at  least  the  software  designer  need  not  be  concerned  about  the  merits 
of  the  underlying  theoretical  procedure. 

Such  finite  decision  procedures  do  not  exist  for  most  optimization 
problems  of  the  form  PI  , and  it  is  necessary  to  make  some  assumptions 
about  the  properties  of  F and  order  to  have  any  hope  of  solving 

these  problems  with  an  automatic  routine.  The  creation  of  an  optimization 
library  thus  involves  not  only  programming  considerations,  but  also  defini- 
tion of  a model  of  the  problems  to  be  solved  and  selection  of  appropriate 
algorithms . 
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I'hc  objective  of  the  present  library  it  to  provide  sound,  careful 


imi)lemcntations  of  reliable  methods  for  solving  useful  categories  of 
optimization  problems.  The  routines  arc  written  in  a portable  subset  of 
ANSI  Fortran  (ANSI,  . so  that  they  may  be  used  with  a minimum  of 

change  on  several  machine  ranges  The  library  is  continually  modified 
to  represent  as  closely  as  possible  the  prevailing  "state  of  the  art"  in 
numerical  optimization. 

2 .2  Choice  of  algorithms 

The  efficient  solution  of  optimization  problems  by  a single,  all- 
purpose, method  is  not  possible,  and  the  structure  of  the  library  reflects 
the  variety  of  theoretical  algorithms.  This  diversity  is  not  surprising, 
in  view  of  the  complexity  of  optimization  problems.  Although  solving  the 
eigenvalue  problem  is  based  in  most  cases  on  a finite  decision  procedure, 
the  EISPACK  set  of  eigensystem  software  (Smith  et  al,,  19711b)  contains  'yPt 
distinct  explicitly  called  subroutines  and  12  driver  subroutines,  the 
choice  of  which  depends  on  properties  of  tlie  given  problem;  the  number  of 
alternative  routines  for  optimization  problems  will  necessarily  be  even 
larger . 

Optimization  algorithms  are  designed  to  solve  particular  categories 
of  problem,  where  each  category  is  defined  by  properties  of  the  objective 
and  constraint  functions,  as  illustrated  below. 
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Properties  of  (c 
No  constraints 
Upper  and  lower  bounds 
Linear 

Sparse  linear 
Nonlinear 

The  choice  of  algorithm  depends  not  only  on  the  type  of  problem,  but 
also  on  the  available  information  about  the  problem  functions,  the  dimen- 
sion of  the  problem,  the  cost  of  evaluating  the  problem  functions,  and  so  t-r . 

The  choice  of  which  algorithms  to  include  in  the  library  is  complicated, 
in  part  because  of  the  difficulty  of  assessing  the  relative  merits  of  dif- 
ferent algorithms  that  solve  the  same  category  of  problems.  Any  comparison 
among  algorithms  inherently  includes  a merit  function  to  be  maximized,  and 
for  optimization  problems  there  is  no  universal  agreement  as  to  what  con- 
siderations should  be  included  in  the  merit  function.  Even  for  relatively 
straightforward  problems  in  linear  algebra,  such  as  solving  a linear  system, 
it  is  possible  to  create  examples  for  which  the  generally  accepted  "best" 
algorithms  are  non-optimal.  For  the  more  difficult  problem  of  nonstiff 
ordinary  differential  equations,  Shampine,  Watts  and  Davenport  (I9TL)  pro- 
vide a lengthy  and  detailed  discussion  of  criteria  for  measuring  the  per- 
formance of  ODE  codes,  and  repeatedly  emphasize  that  the  decision  about 
wliich  code  is  "best”  depends  on  the  problem  to  be  solved  and  on  the  quali- 
ties that  are  most  important  to  the  particular  user.  A similar  statement 
can  be  made  concerning  optimization  problems;  the  determination  of  the  "best" 


al^’Oiithm  lor  most  problems  depends  on  -i  comp  I ieatin!  analysis  of  the  user's 
needs,  as  well  as  on  problem-dependent  ini oimarion,  much  ol  which  nvay  be 
unkno^^n . Accordingly,  the.  pre.sont  optimi  eat  i on  liurary  is  designed  to  in- 
clude implementations  of  the  algorithms  that  have  performed  most  success- 
fully in  extensive  testing  on  problem.s  that  are  believed  to  be  typical  of 
those  encountered  in  a general  environment  i where  "success"  is  measured  by 
some  reasonable  criterion;  see  Gill  and  Murray,  107^-h,  Chapter  9,  for 
addi t ronal  comments } . 

2 .5  Algorithm  standards 

All  theoretical  algorithms  to  be  included  in  the  library  are  required 
to  be  reliable  and  robust  (see  Cody,  1975,  for  definitions  of  these  terms). 
Any  algorithm  to  be  implemented  for  general  use  should  be  able  to  test 
whether  the  assumptions  upon  which  it  is  based  seem  to  be  satisfied,  and  to 
deal  satisfactorily  with  situations  not  satisfying  those  assumptions.  As 
an  illustration  of  this  standard,  a widely  discussed  issue  in  unconstrained 
optimization  has  been  the  procedure  that  should  be  followed  in  a Newton-type 
algorithm  when  the  Hessian  matrix  at  the  current  point  is  not  positive 
definite  'see  Murray,  1972,  Chapter  1,  for  a detailed  discussion).  Any 
Newton-type  routine  in  the  library  must  be  able  to  detect  an  indefinite 
Hessian  and  carry  on  in  a reasonable  and  stable  manner.  Similar  safe- 
guards should  be  included  in  the  definition  of  every  theoretical  algorithm 
proposed  for  the  library. 

An  additional  algorithmic  consideration  is  that  many  alternative  com- 
putational procedures  exist  which,  although  theoretically  equivalent,  dis- 
play widely  differing  numerical  behavior.  The  distance  between  the  state- 
ment of  a mathematical  process  and  its  implementation  in  finite  precision 
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has  been  emphasized  by  several  authors  (see,  for  example,  Smith  et  al,, 
197^ia).  In  the  present  library,  it  is  considered  essential  to  verify  that 
every  procedure  to  be  carried  out  has  satisfactory  numerical  properties. 
For  example,  in  all  library  routines  for  quasi-Newton  methods  for  uncon- 
strained minimization,  the  Hessian  approximation  is  represented  in  terms 
of  its  Cholesky  factorization,  to  allow  strict  control  over  whether  the 
matrix  is  numerically  positive  definite  (Gill  and  Murray,  1972);  this 
facility  does  not  exist  if  an  approximation  to  the  inverse  Hessian  is 
maintained , 

2 .4  Structure  of  implementations 

The  structure  of  the  library  has  been  chosen  to  reflect  as  much  as 

t 

possible  the  ideas  of  good  design  that  are  sometimes  called  "structured 
programming"  or  "step-wise"  algorithm  design  (see  Dijkstra,  I972).  In 
essence,  an  abstract  algorithm  is  decomposed  into  steps,  each  of  which  is 
successively  sub-divided  until  the  nature  of  the  algorithm  is  revealed  by 
a highly  structured  combination  of  concise  and  well-defined  computations. 
After  this  analysis  has  been  carried  out  for  the  theoretical  algorithm, 
the  computer  implementation  is  designed  to  display  the  same  structure. 

In  our  view,  such  a systematic  analysis  of  algorithms  is  essential  to 
obtain  good  numerical  software,  and  should  be  completed  before  a single 
line  of  code  is  written. 

The  process  of  dividing  an  algorithm  into  components  of  decreasing 
complexity  has  been  called  many  names  — for  example,  the  term  "modulariza- 
tion" is  sometimes  applied  in  this  way.  Unfortunately,  "modularization" 
is  often  used  to  imply  an  analysis  of  the  components  of  an  "algorithm" 


that  occurs  alter  a >.  onriuit  cr  ir.ij)  leiiu'utati  on  .ilroiuiy  exists  — the  reverse 
of  the  proceiiure  described  above  (see  Miukott,  I'l'f-j).  ihereCore,  we  wish 
to  avoid  the  tenn  "modal arizat i on",  becaust'  an  after-the-iact  analysis  of 
computer  proijrmiis  is  contrary  lo  the  spirit  ot  ihe  desired  approach.  We 
preler  instead  to  define  the  aU-orithms  and  relat  id  imp  I eionntations  in  the 
present  1 ibrairy  as  "use.iullv  st  met  ur.-,d'  ;ui  jksti.i,  197b,  p.  9)- 

One  of  the  ’leaefits  of  this  bind  of  analysts  is  that  it  is  possible 
to  identity  closely  related  or  identical  sub-steps  in  a "ariety  of  algo- 
rithirus;  hence,  the  same  subroutine  can  be  use.d  in  numerous  distinct  con- 
texts. A further  advantage  is  that  this  approach  n.ay  reveal  structural 
similarities  in  algorithms  that  were  considered  initially  to  be  quite 
different,  ba^auc"  their  underlving  nature  was  disguised  by  notation  or 
computational  detail  (see  the  examples  given  by  Lawson,  I9T'  )• 

Well -structured  programs  are  vital  if  extensions,  modifications,  or 
genera  I illations  are  to  occur  without  undue  complication . Since  active 
rescar'.n  is  continuing  in  all  areas  of  optimization,  it  is  probable  that 
improvements  will  be  made  to  various  aspects  of  every  algorithm.  If  the 
structure  of  the  underlying  algorithms  has  been  accurately  reflected  in 
the  computer  implementations,  it  .should  be  easy  to  extend  the  algorithmic 
capabilities  of  the  library.  For  exatT!;)le,  a new  quasi-Newton  update 
formula  can  be  included  as  a minor  alteration  of  an  existing  quasi-Nev.'ton 
routine . 

Finall . carclul  structuring  facilitates  the  changes  that  might  result 
trom  a n'ore  powerful  progranmdng  language,  or  ' rom  a generalized  data  struc- 
ture provision.  For  example,  if  it  became  possible  to  store  symmetric 
matrices  in  minimal  storage  with  i;fficient  double  indexing,  luany  schemes 
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devised  to  avoid  wasting  storage  under  the  current  limitations  of  Fortran 
(see  Section  "5.?)  would  become  unnecessary;  tlie  relevant  routines  should 


be  designed  so  that  such  a change  would  have  a localized  and  well-understood 
effect.  Good  structuring  can  avoid  the  serious  danger  of  "thinking  in 
Fortran"  rather  than  perceiving  the  overall  nature  of  an  algorithm;  a con- 
cern with  the  details  of  Fortran  often  leads  to  primitive  programming  prac- 
tices, which  then  perpetuate  the  restrictions  of  present  software  tools  and 
make  it  virtually  impossible  to  take  advantage  of  any  improvements. 

2 .5  Structure  vs  program  efficiency 

A further  principle  of  the  library  design  is  that,  if  a choice  is  neces- 
sary, good  structure  rather  than  "efficiency"  will  be  emphasized.  The  term 
"efficiency"  is  taken  here  to  apply  only  to  programming  considerations;  the 
efficiency  of  the  numerical  procedures  to  be  followed  is  included  in  the 
factors  that  determine  the  choice  of  algorithm. 

The  first  point  to  be  stressed  is  that  the  efficiency  of  a software 
design  should  be  viewed  within  the  full  context  of  the  total  amount  of  effort 
required  to  develop  the  routines  as  well  as  solve  the  given  problems,  and  not 
(as  it  often  is)  defined  solely  in  terms  of  one-time  execution  speed.  The 
efficiency  of  a library  design  involves  not  only  the  cost  of  an  isolated 
execution  of  a few  subroutines,  but  also  the  tasks  of  verification,  testing, 
and  modification.  It  is  widely  accepted  that  the  need  to  create  a manageable 
and  reliable  set  of  routines  should  dominate  a concern  with  run-time  efficiency: 

"Premature  emphasis  on  efficiency  is  a big  mistake  i^ich  may  well  be  the 
source  of  most  programming  complexity  and  grief  ....  we  should  strive  most  of 
all  for  a program  that  is  easy  to  understand  and  almost  sure  to  work" 

(Fnuth,  197't,  p.  29)1)- 
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A so.'.oni)  coiisiderat-ion  is  thaL  Llie  proco.ss  o[  writing  "efficient" 
Fortran  coc!n  is  not  at  all  straiglitf  orward . Given  tliat  library  routines 
are  written  in  a portable  subset  of  ANSI  Fortran,  to  be  run  with  a minimum 
ot  change  on  different  computers,  it  is  im])Ossible  in  many  cases  to  write 
code  that  is  effieicnt  on  all  machines.  The  surprising  examples  cited  by 
l'arli!tt  and  Wang  ii/f'')  oi  fei  a decisive  illustration  that  the  effort  to 
write  "efficient''  Fortran  requires  a detailed  a priori  knowledge  of  the 
compile.r  and  nuichine  to  he  used.  In  fact,  someone  vdio  has  made  an  attempt 
to  write  "optimized''  Fortran  may  find  that  tl>e  code  written  by  a less 
sophisticated  programmer  actually  runs  faster. 


3 . Discussion  of  PesiKn  Decisions 
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3 . 1 Communication  by  parameter  List 

A fundamental  early  decision  concerns  the  mechanism  by  which  the  various 
library  subroutines  communicate  information  to  one  another  and  to  the  user. 
Because  the  routines  are  written  in  standard  Fortran,  the  two  possible  means 
of  communication  are:  (1)  through  the  formal  parameter  list;  and  (!•)  through 

COMMON  storage.  Some  previous  optimization  codes  have  chosen  alternative  (2) 
fLasdon  et  al.,  I9TO;  Mylander,  Holmes  and  McCormick,  l97^t ; Rosen  and  Wagner, 
19T5)>  while  others  have  chosen  (1)  (Fletcher,  1972);  however,  the  reasons 
for  the  decision  were  not  specified  in  any  of  these  cases. 

For  the  present  library,  it  was  decided  to  restrict  communication  among 
routines  entirely  to  the  parameter  list  mechanism.  COMMON  storage  is  some- 
times viewed  as  the  Fortran  equivalent  of  global  variables  in  a block- 
structured  language,  but  the  use  of  COMMON  suffers  from  many  defects  not 
present  with  the  facilities  of  true  glooal  variables.  We  shall  give  three 
reasons  for  the  total  omission  of  the  use  of  COMMON  by  library  routines: 

(1)  With  COMMON,  it  is  impossible  to  include  arrays  of  variable 
dimension.  In  most  optimization  applications,  the  lengths  of  the  program 
arrays  depend  conceptually  on  the  dimension  of  the  problem  — for  example, 
the  vector  of  independent  variables  is  represented  as  a floating-point 
vector  of  length  n . Large  fixed-length  arrays  in  COMMON,  of  sufficient 
size  to  solve  most  problems,  would  result  in  a serious  waste  of  storage, 
particularly  since  two-dimensional  arrays  are  often  involved.  Although 
the  user  could  be  asked  to  re-dimension  all  arrays  to  correspond  to  the 
size  of  his  problem,  such  a requirement  makes  an  unreasonable  demand  on 
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the  liSv.r , ^u•ul  also  rules  out  (■.lia  possil'iiity  lor  the  library  rouLlnes  to  be 
pre-cop,ipi  leii . 

(;')  It  has  been  our  eKperiencc  that  otten  the  user  has  already  set  up 
the  COMMON  area  within  his  program  before  h.e  is  ready  to  use  the  library 
routines.  In  such  cases,  the  necessity  to  alter  COMMON  storage  throughout 
the  program  is  e.xtreuiely  inconvenient,  and  frequently  leads  to  obscure  and 
complicated  programming  errors. 

(y)  The  final  justification  for  coninunication  solely  by  parameter 
list,  which  underlies  the  two  previous  reasons,  is  that  this  restriction 
is  a consequence  of  the  ideas  of  good  structure  discussed  in  Section  f.b. 

The  ability  to  understand  a given  subroutine  is  greatly  enhanced  when  all 
the  quantities  upon  which  its  execution  depends  are  in  the  parameter  list. 

In  addition,  the  limitation  to  parameter  list  communication  means  that  the 
internal  logic  of  each  subroutine  tan  be  made  independent  of  the  external 
environment  into  which  it  is  to  be  placed.  Since  all  infornation  is  passed 
through  strictly  controlled  channels,  there  is  no  worry  that  the  subroutine 
will  produce  unwanted  side-effects  'for  examjjie,  by  overwriting  a portion  of 
COMMON  because  of  an  accidental  coincidence  ol  variable  names). 

The  dimensions  of  any  variable-length  array  are  always  included  in  the 
parameter  list  along  with  the  array  name,  fortran  alloxv's  an  array  parameter 
say  "B',  to  be  declared  within  a subroutine  by  "UIMENSION  regardless 

of  the  true  length  of  the  actual  parameter;  but  such  a short  cut  is  not 
allowed  by  some  compilers  — tor  cxair.ple,  the  WATFIV  compiler  (see  Cress  et 
al.,  1970)  — and  it  may  also  disguise  the  nature  of  the  arrays  so  declared. 
All  library  routines  therefore  include  the  correct  diniensions  as  parameters 
in  the  calling  sequence. 
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j 3 .2  Choice  of  data  structures 

I 

i The  careful  selection  of  data  structures  is  an  essential  part  of  the 

' design  of  computer  programs.  A good  high-level  programming  language  allows 

the  programmer  to  define  a data  structure  so  that  its  manipulation  within 
the  program  resembles  its  original  role  in  the  theoretical  algorithm. 

Ideally,  the  programmer  should  not  be  required  to  transform  convenient 
conceptual  data  structures  to  an  inconvenient  machine-based  representation; 
but  the  highly  limited  capabilities  of  Fortran  necessarily  cause  a depar- 
ture from  this  ideal. 

In  numerical  applications  such  as  optimization,  vectors  and  matrices 
are  the  primary  conceptual  data  structures;  it  might  seem  straightforward 
to  transform  these  into  one-  and  two-dimensional  Fortran  arrays,  respec- 
tively. Such  a transformation  is  convenient  in  some  cases  — for  example, 
the  vector  of  independent  variables  can  be  represented  in  an  obvious  way  as 
a one-dimensional  floating-point  array.  In  other  instances,  however,  the 
choice  may  be  complicated  by  programming  considerations. 

Some  Fortran  programmers  represent  matrices  as  one-dimensional  arrays, 
because  true  double  indexing  does  not  exist  in  Fortran.  The  element  A(I,J) 
is  obtained  as  element  [ (j-1)  * (declared  row  dimension  of  A)  + l]  of  the 
block  of  contiguous  storage  locations  representing  the  matrix  A,  and  thus 
there  would  seem  to  be  an  inherent  overhead  cost  associated  with  accessing 
two-dimensional  arrays.  The  additional  cost  of  processing  a two-dimensional 
array  compared  to  a one-dimensional  array  depends  on  several  complicated  j 

factors,  which  include  the  hardware  features,  the  compiler,  the  array  dimen- 
sions, the  loop  structure  of  the  code,  and  the  particular  indexing  strategy 
chosen . 
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To  obtain  some  estiriLates  of  the  rolativo.  execution  limes  for  theoreti- 
cally equivalent  manipulations  of  one-  and  two-dimensional  arrays,  numerous 
experiments  were,  run  on  an  IBM  -.'70/K'?  computer,  with  the;  Fortran  U compiler, 
opt  - 0 (unoptinized ) and  opt  (optimized).  As  expected,  the  largest 

increase  in  running  time  was  displayed  for  an  unoptimized  double  do-loop  ' 

that  initializes  ali  elements  of  a tv.'O-diraeiisional  array  to  a constant, 
since  the  loop  and  Inde:^  processing  then  dominate  the  trivial  calculation 
within  the  loop;  the  worst  case  i,a  increase  over  the  one-dimensional 

loop  ) occurred  when  the  index  of  the  outer  loop  was  50  times  larger  than 
for  the  inner,  because  of  the  repeated  overhead  in  setting  up  the  inner 
loop.  With  the  opt  -r  0 compiler,  there  were  a few  other  instances  when 
the  combination  of  unoptimized  code  and  an  unfortunate  choice  of  loop 
structure  led  to  a noticeable  ( IC  - 15“^)  increase  in  execution  time  for  a 
particular  section  of  code  involving  Lcie  two-dimensional  array.  However, 
in  all  remaining  cases,  which  involved  a selection  of  typical  array  calcu- 
lations, the  difference  in  execution  time  for  processing  the  two-dimensional 
array  was  negligible  (often  less  than  the  fluctuation  in  timing  estimates), 
especially  with  the  optimizing  compiler.  This  result  confirms  the  inability 
to  predict  the  relative  efficiencies  of  sections  of  Fortran  code  without 
knowledge  of  the  compiler  and  the  particular  data  to  be  processed. 

In  view  of  this  uncertainty,  it  seems  sensible  to  use  the  most  straight- 
forr.'/ard  data  structures,  especially  since  the  representation  of  a m.atrix  as 
a one-dimensional  array  often  causes  confusion  and  complexity,  both  for  the 
user  and  the  library  programmers.  In  the  library  least -squares  routines, 
for  example,  if  the  user  were  required  to  store  the  cojnponents  of  the 
Jacobian  matrix  in  a one-dimensional  array  according  to  some  non-triviul 
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indexing  strategy,  there  would  be  a high  probability  thar  tlie  inionnation 
would  be  jumbled.  The  discarding  of  double  indexing,  is  in  most  cases  an 
excessive  concession  to  the  limitations  oi  current  I'urtran  in  quest  of  the 
dubious  possibility  of  marginally  faster  execution  speed. 

There  are  instances,  however,  when  a theoretical  two-dimensional 
structure  is  represented  in  the  library  routines  by  one-dimensional  arrays.  i 

A frequent  question  in  implementations  of  linear  algebraic  techniques  con-  | 

cerns  a suitable  program  representation  of  a symmetric  matrix  (similar  con- 
siderations also  apply  for  a triangular  matrix).  Conceptually,  a symmetric 
matrix  has  a two-dimensional  character;  but  if  it  is  represented  by  a two- 
dimensional  array,  almost  half  the  storage  is  "wasted"  in  the  sense  that 
duplicated  information  is  stored  that  could  be  found  by  a simple  interchange 
of  indices. 

If  unlimited  storage  were  available,  it  would  be  feasible  to  retain 
the  natural  two-dimensional  representation  for  a symmetric  matrix.  In  most 
applications,  however,  it  is  unreasonable  to  use  0 storage  locations 

to  contain  easily  accessible  duplicated  information,  and  thus  symmetric  I 

matrices  are  often  represented  by  a combination  of  one-dimensional  arrays.  | 

For  example,  in  the  modified  Newton  routines  for  unconstrained  optimization,  I 

the  symmetric  Hessian  matrix  is  stored  in  two  one-dimensicnal  arrays  — the 
diagonal  elements  in  a vector  of  length  n , and  the  off-diagonal  elements 
in  a vector  of  length  nfn-l)/p  (stored  by  rows).  This  representation  is 
extremely  convenient,  because  the  (possibly)  modified  Hessian  matrix  is 
factorized  into  the  form 

- T 

G = L D L , 
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where  D is  a diagonal  matrix,  auu  I, 


is  a unit  lower  triangular  matrix. 


Hence,  the  vectors  that  represent  the  original  matrix  can  be  overwritten 
with  the  vectors  that  represent  the  Cholesky  factors.  The  resulting  data 
structure  is  not  entirely  straightiorvjard , since  a non-trivial  transforma- 
tion of  indices  is  necessary  to  obtain  the  general  off-diagonal  element 


[Gi.)  or 
this  data 
nif icant 
c laritv . 


(L  .]  from  the  one-dimensional  array. 
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structure  is  confined  to  a few  routines, 
savings  in  storage  is  considered  to  ouLwei 


However,  the  effect  of 
and  the  benefit  of  sig- 
gh  the  localized  loss  of 


3.3  Local  arrays 

Optimization  algorithms  typically  involve  numerous  temporary  arrays 
x^hose  dimensions  are  related  to  the  dimension  of  the  problem  to  be  solved. 
The  issue  to  be  considered  in  this  section  is  the  mechanism  by  which  storage 
is  allocated  to  such  local  arrays.  In  Algol,  a procedure  may  contain  decla- 
rations of  variab le- length  local  arrays,  where  the  dimensions  depend  on  pa- 
rameters of  the  procedure;  but  in  Fortran,  locally  declared  arrays  must  be 
of  fixed  dimension.  The  choices  for  providing  temporary  arrays  in  Fortran 
subroutines  thus  arc: 

(1)  declare  the  arrays  locally  to  each  subroutine,  with  fixed  dimen- 
sions large  enough  tor  most  anticipated  problems.  The  user  must  modify  the 
source  text  of  every  such  routine  if  the  declared  dimensions  are  insuf- 
ficient ; 

f:  ) declare  the  space  for  any  local  arrays  outside  the.  subroutine, 
rh'  vectors  that  contain  this  space  are  then  included  in  the  parameter  list 
of  the  subroutines  that  require  te.mporarv  arrays,  along  with  the  appropriate 
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intagars  Chat  dcEino.  the  dimensions  (see  Section  .1). 

Alternative.  H, } has  the  same  disadvantages  discussed  in  Section  i.l 
with  respect  to  the  use  of  COMMON  storage,  'ihere  is  a significant  waste 
of  storage  if  large  arrays  of  fixed  di-mansion  are  declared  in  many  sub- 
routines^ especially  since  the  locally  declared  storage  cannot  be  shared 
by  several  subroutines;  and  it  is  inconvenient  to  require  the  user  to 
modify  local  array  declarations  tiiroughout  the  library  routines.  For 
some  types  of  numerical  problems,  the  waste  of  storagi;  due  to  fixed-size 
locally  declared  arrays  may  not  be  serious  enough  to  justify  consideration 
of  alternatives.  For  example,  in  the  ordinary  differential  equation  code 
DE  (Shampine  and  Gordon,  1975))  fixed  local  arrays  are  acceptable  because 
the  dimensions  of  the  arrays  are  known  in  general  to  be  of  reasonable  size, 
and  the  maximum  amount  of  unused  storage  is  relatively  insignificant  com- 
pared to  the  program  size  (since  all  arrays  except  one  are  one-dimensional). 
In  optimization  routines,  however,  large  two-dimensional  arrays  are  fre- 
quently required,  and  the  problem  of  wasted  storage  cannot  be  ignored. 
Consequently,  alternative  (2)  i.s  used  in  the  present  library. 

The  decision  to  avoid  declaring  local  arrays  whose  conceptual  dimension 
is  problem-dependent  introduces  into  the  computer  implementations  a deviation 
from  the  straightforward  framework  of  the  original  algorithms,  since  the 
parameter  list  of  a given  subroutine  will  include  not  only  quantities  that 
represent  its  input  and  output,  but  also  quantities  whose  role  is  purely  to 
provide  storage  for  temporary  intermediate  results.  This  disadvantage  is 
nontheless  believed  to  be  outweighed  by  the  increase  in  flexibility  and  the 


savings  of  storage. 


Libfiiry  subrout  iiieb  hicl  udu  il  L var iub  Lu-s i ze  temporary  arrays  and 
their  diuiensions  in  ttie  calling  sequence,  i see  Section  li  such  local 

arrays  are  declared  individually,  the  f'L'j,<,Lh  in  the  length  oi  some  parameter 


lists  can  become  extreme,  ihere  is  .10  inherent  objection  t.o  long  paramo.ter 
lists,  il  each  parameter  has  a signiiicant  role  in  the  definition  of  the 
input  or  output  of  the  subroutine.  However,  this  qual il ication  docs  not 
i.mply  that  there  must  .lecessarily  be  a ; eparate  parameter  name  for  every 
variable- length  temporary  array,  whii  li  has  a pux'ely  internal  significance. 
The  EISPACK  subroutine  TSTURM  has  six  temporary  vectors  in  its  calling 
sequence,  and  this  number  may  not  seem  excessive;  but  some  optimization 
routines  have  well  over  hO  distinct  temporary  arrays.  The  clutter  of  a 
long  sequen.,e  of  temporary  array  names  will  tend  to  distract  from  the 
insight  into  the  nature  of  a subroutine  that  the  elements  of  the  formal 
parameter  list  should  give. 

An  alternative  to  including  individual  work-space  parameters  in  the 
calling  sequences  is  provided  by  partitioning  larger  ai'rays . It  is 
acceptable  in  standard  Fortran  to  call  a subroutine  'v'hose  k-th  formal 
parameter  is  an  array  name  with  an  element  of  an  array  as  the  k-th  actual 
parameter.  This  device  is  consistauL  because  an  array  is  completely  de- 
fined in  Fortran  by  tlie  address  of  its  iiist  element  and  an  indexing  rule. 
Thus,  an  array  parameter  is  defined  within  a subroutine  as  the  set  of 
contiguous  memory  locations  whose  first  element  is  at  the  address  specified 
by  the  corresponding  actual  parameter.  This  convention  lor  specifying  an 
.irriy  means  tliat  one  large,  array  can  be  p.irt it ioriud  into  several  smaller 
sub-umys.  The  location  oi  the  first  element  01  each  partition  cun  then 
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serve  in  a subroutine  caii  a.;  die  iutual  lAuameter  eni-respiauling  to  an 
array  of  appropriate  iengtli,  since  the  execution  of  the  subroutine  gen- 
erates accesses  only  to  eiemants  of  the  Jesirc^l  sub-vector. 

The  flexibility  provided  by  work  si>ace  arrays  can  be  exploited  if  the 
library  routines  are  developed  and  exteitded.  lor  example,  ii  a new  method 
is  devised  for  solving  a givt- i opti mixal i on  prohlei.i,  which  requires  more 
intermediate  calculations  than  the  original,  the  form  of  the  existing 
library  module  can  be  retained  while  implementing  the  new  algoritiim,  by 
simply  using  additional  elements  of  the  already  ]>resent  work  space  arrays. 

Tlie  technique  of  partitioning  one  array  into  many  can  lead  to  a lack 
of  clarity  in  the  resulting  code,  since  the  names  of  the  conceptual  arrays 
of  the  original  algorithm  do  not  appear.  'flie  library  routines  apply  this 
solution  only  with  caution,  and  in  every  instance  the  code  is  surrounded 
by  detailed  comments  explaining  the  partitioning. 

Tlie  difficulty  may  eventually  be  resolved  by  the  use  of  pre- processors 
that  make  intelligent  text  substitutions.  With  such  a pre- processor,  all 
textual  references  to  the  desired  local  array  name  could  be  replaced  before 
compilation  by  references  to  suitable  elements  of  the  work  space  array, 
i.e.,  if  the  vector  P , of  length  N , began  in  location  Nil  of  the 
vector  W , references  to  could  be  replaced  by  W( N i I)".  .Mtcrna- 

tively,  the  recent  ''MA.P"  proposal  to  modify  the  ANSI  Fortran  standards 
(IFIP  Working  Group  cl.'j,  197'0  provides  a similar  means  of  gaining  such 
flexibility. 


■j.  j Intel  £acc  wiili  tlie  ur.oi: 

In  dt-^igiu'ng  library  software,  it  niisi  be  iii-ci')eJ  e.'iiliiilly  liow  much 
uC  the  llu-.ory  and  ciOTp  1 exi ty  uf  a nu’iericai  algoridiin  t/ie  user  need.s  to 
understand  beiore  making  use  of  the  associ.ited  routines.  liie  controversial 
statomeat  that  a possible  purjKise  ot  mathematical  sofiware  is  to  'relieve 
[the  user]...  of  any  need  to  think"  ( Davis  and  Rabinowitz, 
in  the  context  of  automatic  quadialurc'  is  confirmed  as  tlie  ideal  of  many 
users  by  other  authors  (see  Dixon,  1?7'' I Kahaner,  iui'l;  Kuki,  19Y1),  and 
by  the  experience  of  an.-one  who  has  ever  acted  as  a consultant  concerning 
numerical  problems.  Because  many  users  found  it  burdensome  to  have  to  call 
several  (even  tv/o'l  subroutines  that  represent  the  conceptual  steps  of  an 
eigensys ceiii  computation,  the  LISPACK  system  now  provides  driver  routines 
which  call  all  the  related  subroutines  required  for  a particular  problem. 
Many  people  refuse  to  use  routines  that  require  more  than  a marginal  effort 
to  understand,  and  wish  to  know  nothing  of  the  possible  subtleties  involved 
in  solving  tlie  problem.  A consequence  of  this  attitude  among  users  is  that 
many  software  designers  try  to  keep  all  calling  sequences  as  short  as  pos- 
sible (see  Kahaner,  I'YTl;  Bhampine  and  Gordon,  because  users  are 

reluctant  to  try  subroutines  witii  long  calling  sequences  and  may  even  prefer 
an  inferior  routine  with  a shorter  parameter  list.  Hiis  situation  leads  to 
the  feeling  that  software  designers  shc<uld  try  tu  remove  any  hint  of  dif- 
ficulty for  the  user. 

i\n  opposing  view  is  that  tlic  user  should  be  forced  to  become  at  least 
minimally  educated  and  informcel  about  the  numerical  procedure.s  to  be  carried 
out  by  a library  routine.  It  can  be  argued  that,  by  catering  to  tlu:  user's 


desire  for  simplicity,  the  authors  of  matliemat  icaJ  s<j)  tware  are  performing 
a disservice  if  the  form  of  the  software  gives  the  false  Jiiipression  that  a 
problem  is  easy  and  ^;trai;’htforward  to  solve.  In  particular,  an  atcemjjt 
to  enforce  short  parajiieter  lists  to  please  users  may  compromise  the  per- 
formance of  an  optimi;^at  ion  algorithm,  because  the  adiustment  of  selected 
tolerances  and  convergence  criteria  can  have  a significant  influence  on  the 
efficiency,  and  even  the  success,  of  a given  metliod  applied  to  an  optimiza- 
tion problem.  'I'he  most  suitable  values  for  these  parameters  typically 
depend  on  knowledge  that  only  the  user  has,  which  would  seem  to  imply  that 
for  optimum  performance  the  library  routines  should  contain  all  such 
parameters  in  the  calling  sequence. 

The  dilemma  in  library  uesign  is  that  the  creation  of  superb  routines 
will  be  unproductive  if  no  one  uses  them  because  of  their  excessive  com- 
plexity; on  the  other  hand,  overzealous  simplification  to  satisfy  some 
users  may  lead  to  the  loss  of  highly  desirable  flexibility  for  others. 

Rather  than  compromise  between  flexibility  and  ease  of  use,  the  present 
library  has  been  designed  to  satisfy  both  objectives  by  providing  two 
routines  — a general  and  a "simple''  — for  each  purpose.  A similar  idea  is 
seen  in  the  EISPACK  collection  (Smith  et  al.,  I97^b)  and  in  the  ODE  rou- 
tines described  by  Shampine  and  Gordon  (19'/'5)-  If  a particular  parameter 
can  be  chosen  automatically  (given  that  certain  assumptions  are  satisfied), 
the  uncertain  or  uninterested  user  may  choose  a "simple"  or  "easy-to-use  " 
routine,  for  which  only  a minimum  number  of  parameters  need  to  be  specified. 
However,  the  possibility  is  retained  that  a well-informed  user  can  enhance 
the  performance  of  an  algorithm  by  exploiting  his  knowlcilge  of  the  problem 


to  select  tho  appropriate  paraincters  ol  tlie  general  routine.  'AH  user.s  are 
not  required  to  think,  but  tlie  knowledgeable  user  is  at  least  given  the  oppor- 
tunity of  Lliinking.) 


Althougli  this  scheme  would  appear  to  double  tlie  amount  of  code,  the 
easy-to-use  routines  are  simply  "front  end.s",  wliich  call  the  more  general 
routines;  thus,  the  increase  in  code  is  not  significant.  Similarly,  the 
amount  of  documentation  is  not  doubled,  .since  one  of  the  objectives  in  pro- 
viding easy-to-use  programs  is  to  reduce  the  amount  of  information  that  tlie 
user  must  assimilate  before  he  can  solve  his  problem. 

It  should  be  noted  that  a development  in  programming  language  control 
structures  may  soon  make  the  concern  about  long  calling  sequences  unneces- 
sary. This  development  is  the  "keyword"  parameter  communication  mechanism, 
described  in  Hardgrave  f 19'f’C)  and  Zahn  ( 1975)-  The  idea  is  that  the  actual 
parameters  in  a subroutine  call  are  not  specified  by  position,  but  rather  by 
keyword.  The  user  thus  needs  to  include  in  the  actual  calling  sequence  only 
the  parameters  that  he  wishes  to  specify,  and  the  others  are  set  by  default, 
so  that  a flexible  routine  can  also  serve  as  "easy-to-use".  in  fact,  the 
keyv;ord  parameter  mechanism  would  allow  a range  of  alternatives,  from  speci- 
fying no  optional  parameters  to  a full  calling  sequence. 

9 . 6 Service  routines 

The  "simple  ' routines  described  in  Section  9-5  have  been  provided  in 
response  to  user  demand  for  short  calling  sequences.  However,  the  remaining 
parameters,  which  arc;  masked  from  the  user  of  the  simple  routine,  must  none- 
fheies;;  be  chosen.  lliere  are  two  way.s  t<5  automate  the  process  of  parameter 
se  Lee  tion. 


First,  a reasoiiablu  a priuri  clioicc  lor  some  paiamtLers  is  possible  il 
tlie  problem  to  be  solved  satisfies  a specified  set  of  assompl ions . For 
example,  the  clu  ice  of  teriioiiatiou  critv-tra  depends  on  tin  [.robleiid  r.  staling; 
if  the  problem  is  assumed  lo  be  'uell  scaled'  in  a given  sense,  then  Llie 
convorgeiice  criteria  can  bt  based  on  t.lie  machine  (irecision.  Similarly,  a 
"sensible''  value  of  the  parameter  that  controls  terii.inatiim  cf  the  step* 
length  algoritlim  can  be  selected  based  on  wide  computational  experience. 

The  parameters  chosen  in  this  way  are  considered  to  be  'safe",  in  the  sense 
that  the  assumed  properties  ol  the  model  [jroblcm  upon  which  their  selection 
is  based  will  correspond  fairly  well  to  those  of  most  pioblems  to  he  solved. 

However,  certain  other  parameters  are  considered  to  be  too  sensitive 
to  the  given  problem  to  allow  an  a priori  choice.  Consi-Liuen tly,  a second, 
more  complicated,  service  to  the  user  is  provided  by  library  routines  that 
automate  problem- dependent  decisions.  llie  procedures  for  parameter  selec- 
tion that  depend  on  properties  of  the  specific  problem  functions  can  usually 
be  automated  to  some  extent,  and  a carefully  designed  automatic  routine  should 
give  a better  result  than  the  random  guess  chat  might  well  come  from  an  un- 
informed user.  For  example,  to  solve  an  unconstrained  opt iraization  problem 
with  a quasi-Newton  method  based  on  finite-difference  approximations  to  gra- 
dients, a finite-difference  interval  must  be  provided  for  each  variable. 

The  values  of  the  difference  intervals  are  e.flen  critical  in  tl.e  performance 
of  the  algorithm,  but  many  users  have  no  idea  how  to  choose  this  set  oi 
intervals.  It  has  therefore,  been  considered  worthwhile  to  iirovide  a library 
routine  to  determine  an  initial  sensible  set  of  i initc-di  f ferenco  interval.s; 
this  routine  (FRMDEL;  sec  Gill  et  al.,  l^'|-'^c)  calls  the  user- provided  func- 
tion, and  chooses  the  intervals  based  on  examining  the  magniiudcs  oi 


■ aacL' 1 laL  ion  and  irniicalion  error'i  in  Liu-  dirivat  i v<-  approximation.  There 
is  no  I'liarantDod  tt'chniq.ie  for  dif  fertpc*;- interva  I selection,  because  the 
values  obtained  at  one  point  may  be  unrc  i ,o:iab  1 e at  a different  point  if 
the  pro!)  iom- func  t ion ' s bidiaylor  alter.s  drastically.  if,  hov/ever,  the 
scaling  of  the  problem  does  not  change  too  much  as  the  minimization  proceeds 
(a  situation  that  usually  holds  in  practice),  the-  initial  .set  of  intervals 
remains  a gooel  choice.  bucli  "service  ' routiners  in  the  library  represent 
tile  view  that  it  is  better  to  automate  parar.ieter  selection  in  a sensible 
way  if  the  user  cannot  provide  the  infornut 1 on,  rathee  than  incur  the  risk 
of  significantly  reduced  efficiency. 

Other  library  subroutines  provide  a further  service  to  the  user  — a 
consistency  check.  One  of  the  most  common  difficulties  in  use  of  opti- 
mization routines  is  that  the  user's  subroutines  incorrectly  evaluate  the 
relevant  partial  derivatives.  Because  exact  gradient  information  normally 
enhances  efficiency  in  all  areas  of  optimization,  the  user  should  be  en- 
couraged to  provide  analytic  derivatives  whenever  possible.  However, 
mistakes  in  the  computation  of  derivatives  can  result  in  serious  and 
obscure  run-time  errors,  as  well  as  complaints  that  the  library  routines 
are  incorrect.  Consequently,  library  routines  are  provided  to  perform  an 
elementary  check  on  the  user-supplied  gradients.  Because  the  checking 
procedure  will  not  be  cost-effective  if  it  requires  too  many  evaluations 
of  user- supplied  functions,  the  present  library  contains  subroutines  that 
compute  finite-difference  approximations  to  two  projections  of  first  and/or 
second  derivatives  (it  would  usually  be  considered  too  expensive  to  compute 
finite  differences  along  all  the  coordinate  directions).  Tlie  approximated 
values  are  compared  with  the  supposedly  exact  values,  and  should  display 
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reasonable  agreemenl.  If  they  do  not,  then  either  the  user  has  made  an 
error  or  the  problem  is  very  badly  scaled  ~ in  both  cases,  the  user  should 
'.ike  corrective  action.  Of  course,  as  with  the  routine  to  select  the 
finite-difference  intervals,  such  a procedure  can  not  guarantee  correct- 
ness; even  if  the  quantities  display  some  agreement,  it  may  simply  be  due 
to  a fortuitous  starting  point.  However,  the  simple  process  described 
above  detects  the  vast  majority  of  user  errors  in  computing  partial  deriva- 
tives. 

5. Y Auxiliary  routines 

i An  "auxiliary  routine"  is  a subroutine  or  function  that  is  not  called 

I 

directly  by  the  user,  but  only  by  other  library  routines, 

1 Subroutines  vs  in-line  code 

The  process  of  analysis  mentioned  in  Section  2.h  leads  to  the  descrip- 
tion of  an  algorithm  as  a structured  and  logically  coherent  collection  of 
computations,  which  must  be  transformed  into  Fortran  code.  Hie  topic  to 
be  discussed  here  is  the  process  of  deciding  whether  a particular  computa- 
tion or  set  of  computations  should  be  implemented  in  the  library  as  a 
separate  subroutine,  or  as  part  of  the  code  body  of  another  routine. 

Some  people  favor  in-line  code  as  a general  practice,  based  on  the 
argument  that  efficiency  (in  the  restricted  sense  of  run-time  speed)  is 
degraded  by  the  overhead  cost  associated  with  a large  number  of  sub- 
routine calls.  In  order  to  achieve  maximum  clarity,  however,  the  division 
of  the  library  into  subroutines  would  be  based  directly  on  the  structure  of 
the  theoretical  algorithms,  with  a one-to-one  relationship  between  the  con- 
ceptual steps  of  an  algorithm  and  those  of  the  corresponding  implementation. 


This  latter  policy  does  not  significantly  decrease  efficiency  in  most  cases 
because  the  execution  time  for  subroutines  tends  to  be  dominated  by  the 
actual  computations  (arithmetic  operations  and  memory  accesses),  and  the 
time  required  to  set  up  the  calling  mechanism  is  negligible  by  comparison. 
However,  in  some  instances  there  is  a marginal  balance  between  the  struc- 
tural benefits  of  creating  an  auxiliary  routine  and  the  consequent  loss  of 
run-time  speed. 

The  following  experiments  were  designed  to  investigate  the  execution 
time  associated  with  a subroutine  call.  First,  code  was  written  to  compute 
in-line  the  scalar  product  of  the  vectors  A and  B,  of  length  N , and  assign 
the  value  to  a variable.  A subroutine  DOTl  (N,  A,  B,  DPROD)  was  written 
to  perform  the  same  calculation,  and  assign  the  result  to  the  variable 
DPROD.  On  an  IBM  370/168,  with  the  Fortran  H compiler,  opt  = 0 and 
opt  = 2 , the  in-line  code  and  calls  to  DOTl  were  executed  3OOO,  toOO,  and 
5000  times,  for  N = 1,2, 10,50, 75> 100,  with  each  run  repeated  twice  to 
observe  the  fluctuation  in  timing  estimates. 

Table  1 gives  the  average  ratio  of  the  execution  times  of  DOTl  to 
those  of  the  in-line  code. 
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1.00 
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3.61 
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Second,  in-line  code  and  a subroutine  l)0r<2(  N,A,  B,  DPR0D,X1,X1;, , . . ,X12) 
were  written,  to  compute  the  scalar  product  of  the  vectors  A and  B as 
before,  and  also  to  assign  the  values  of  some  simple  arithmetic  expressions 
to  the  scalar  variables  X1,X2, . . . ,X12.  llie  in-line  code  and  calls  to  D0T2 
were  executed  repeatedly,  exactly  as  in  the  previous  case.  Table  2 gives 
the  average  ratio  of  execution  times  for  D0T2  to  those  of  the  equivalent 
in-line  code. 


N 

1 

2 

10 

50 

75 

100 

1.89 

1.78 

1.03 

1,00 

1.00 

2. 10 

2.04 

1. 68 

1.21 

1. 14 

1. 10 

Table  2 


It  should  be  stressed  that  the  values  in  Tables  1 and  2 were  computed 
in  a multi- programming  environment,  and  are  inevitably  subject  to  a lack 
of  high  precision;  however,  the  ratios  for  the  various  cases  fluctuated 
only  in  the  hundredths  place,  so  that  they  are  reasonably  accurate  for 
this  machine.  One  must  always  be  cautious  in  drawing  general  conclusions 
from  the  results  of  a particular  timing  experiment,  which  depend  on  the 
compiler  as  well  as  the  machine.  Nonetheless,  the  results  correspond  to 
a priori  expectations,  and  allow  at  least  a rough  measure  of  the  desired 
phenomenon. 

For  these  examples,  the  increase  in  execution  time  attributable  to  a 
subroutine  call  is  very  substantial  for  small  N , noticeable  for  moderate  N, 
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and  ttssontiaily  nef,  1 iy.ib  la  for  lar^o  . ilia  signil  j aance  of  these  dif- 
feraiicas  depends  on  Lha  roia  uL  liia  >^ivtn  cjlculation  within  the  library. 

If  the  computations  of  Ll.a  autira  Library  wari  primarily  a sequence  of 
innar  products  for  small  or  moderate  vaJaa,-i  of  N , it  would  probably  not 
be  worthwhile  to  create  an  auxiliary  dot-|)roduct  subroutine,  since  the 
repeated  subroutine  calls  w-oul.l  clear’y  have  a serious  impact  on  execution 
spaed.  Howe,ver,  in  nonlinear  optimization  most  subr<->utines  perform  order 
I m‘')  or  order  ; h'  ) computations,  and  tb.e  user- supp  1 ied  routines  are 
typically  even  more  time-consuming.  A scalar  product  represents  the 
sin.plost  calculation  that  riiight  be  isolated  into  a subroutine  within  a 
general  optimization  library,  so  that  these  experiments  indicate  the  worst 
degradation  in  execution  tinic  that  could  possibly  occur  as  the  result  of  a 
decision  to  create  an  auxiliary  routine.  Consequently,  the  principle  stated 
earlier  of  dividing  into  subroutines  whenever  possible  to  reflect  the  under- 
lying algorithmic  structure  should  cause  a negligible  loss  of  execution 
speed,  and  a great  gain  in  clarity. 

" . ' . C'  Fle.xibility  of  auxiliary  routines 

It  is  typical  of  optimization  algorithms  that  identical,  or  closely 
related,  computations  must  be  carried  out  in  metliods  to  solve  different 
problems,  in  alternative  methods  to  solve  the  same  problem,  and  in  dif- 
ferent steps  of  a single  method.  If  in-line  code,  or  different  sub- 
routines, for  similar  computations  can  be  replaced  by  a call  to  a single, 
flexible  auxiliary  routine,  tlic  amount  of  source  code  in  Che  library  can 


be  reduced. 


Tlia  following  difficulties  must  be  considered  in  deciding  whether  to 
create  a single  auxiliary  routine  to  carry  out  related  but  slightly  dif- 
ferent computations: 

; a)  the  length  of  the  calling  sequence  may  increase; 

vb)  the  source  code  body  of  the  routine  may  be  larger; 

I c)  the  execution  time  of  the  routine  may  increase  because  of  addi- 
tional logic  or  computation  to  provide  flexibility  (e.g.,  testing 
of  parameter  values). 

Concerning  (a),  the  following  experiment,  similar  to  that  described 
in  Section  was  designed  to  obtain  some  measure  of  the  increase  in 

execution  time  due  to  longer  calling  sequences.  A subroutine 
DOT-’(  N, A,  B,  DPR0D,Xl,X2, . . . ,X1L") , with  16  parameters,  was  written  to  compute 
the  scalar  product  of  the  vectors  A and  B , and  to  assign  this  value  to 
the  variable  DPROD.  The  scalar  variables  Xl, X2, . . . , X12  are  unaltered  by 
DOT;';,  and  thus  the  ratio  of  execution  time  for  DOTJ  to  that  of  DOTl  (Section 
1)  should  indicate  the  effect  on  execution  time  of  12  additional  formal 
parameters. 

Exactly  as  for  the  other  experiments,  the  calls  to  DOTl  and  DOT5  were 
executed  repeatedly  for  various  values  of  N , Table  Z gives  the  average 
ratio  of  execution  times  for  calls  to  DOT;;  to  those  of  DOll. 
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A comi/ar isoii  of  Lliosii  re;julu  with  ’iabtoo  1 and  indicates  that 
tlie  increase  in  cxecntion  t isic  associated  with  additiona!  parameters  in 
the  calliiv'  sequoiue  is  K;e-.  ,:i  .ni  i 1 1 . ai:  t ihan  the  ui',  raise  due  to  culling 
a very  simple  subroutine  in;-tLad  oi  performing  the  calculations  in-line 
(of  course,  the  cautions  nieniioned  iu  Section  ; . | . 1 aliout  timing  tests 
apply  here  as  v/ell).  llic  oveiall  conclusion  is  tliat  long  calling  se- 
quences for  .ii.>;i  1 i arv  routines  will  generally  have  a negligible  effect 
on  execution  speed,  and  the  diificulty  suggested  by  (a)  ’s  not  serious. 

An  analysis  of  di I f i cu I t its  (b)  and  ' c)  will  depend  on  tbe  particu- 
lar case  under  consideration.  In  the  present  library,  flexible  auxiliary 
routines  are  preferreu  to  in-line  code  or  several  similar  subroutines 
because  of  the  simplification  of  library  maintenance,  providing  that  the 
following  conditions  are  satisfied: 

(1)  no  increase  in  the  number  oi  calls  to  user-supplied  routines 
results  from  the  added  flexibility; 

(2)  the  increase  in  the  source  code  body  of  the  routine  is  less 
than  ('approximately)  lO''^; 

( y)  the  increase  in  execution  time  because  of  the  cictra  flexibility 

is  at  least  an  order  of  magnitude  smaller  dian  the  total  execu- 

tion  time  for  the.  routia.i.  For  example,  if  0(  n^ ) operations 

o 

are  performed  in  the  routine,  then  Q( n~)  operations  would  be 
tolerated  to  increase  llexibility. 

■1 . 3 Communication  with  user-supp lied  routines 


For  optimization  problems  of  the  form  Pi  (Section  . .1',  the  user 
must  supply  code  to  evaluate  the  obiective  function,  the  constraint  functions 


and  (possibly)  the  appropriate  derivatives.  Therefore,  the  library  routines 
need  to  contain  a mechanism  for  obtaining  from  the  user  tlie  information 
required  to  proceed  with  the  computation,  and  it  must  be  considered  how 
this  communication  should  take  place. 

An  initial  problem  is  the  actual  control  structure  that  connects  the 
routines  of  the  library  and  che  user.  In  most  similar  applications  (such 
as  quadrature  or  ordinary  differential  equations),  the  user  supplies  sub- 
routines that  define  the  problem  functions,  and  these  subroutines  are 
called  by  the  library  routines;  this  control  structure  is  used  in  the 
present  library.  An  alternative  approach,  used  by  Lawson  and  Krogh  (see 
Krogh,  19^9)5  is  termed  "reverse  communication",  where  the  library  routine 
returns  control  to  the  user  whenever  information  is  required.  A closely 
related  strategy  is  included  in  one  aspect  of  the  library  design  (see 
Section  4.)+),  but  not  at  the  outermost  level. 

The  next  question  that  arises  is  at  what  level  the  library  routines 
should  communicate  with  the  user.  Because  the  library  subroutines  will 
call  the  user- supplied  routines,  the  latter  must  be  defined  with  a fixed 
calling  sequence.  The  composition  of  this  calling  sequence,  and  the  means 
of  interface  with  the  library  routines,  pose  some  interesting  issues,  which 
are  best  illustrated  by  a detailed  example. 

Consider  the  case  of  the  library  subroutine  QNMDER  (Gill  et  al.,  19'f5a), 
an  implementation  of  a revised  quasi-Newton  method  for  unconstrained  opti- 
mization that  uses  exact  gradients  (see  Gill  and  Murray,  1972).  Hie  parame- 
ter list  of  QNMDER  contains  a formal  parameter  SPUN,  a subroutine  which  is 
called  within  QNMDER  whenever  information  concerning  the  function  to  be 
minimized  or  its  gradient  is  required.  It  was  decided  that  the  same 
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bubi'i'u  Li  lu  slu.'iilii  i.iKi.lit.  b.'it.  :he  luncLi'-i  .:iul  i.a<i  i ii , since;  thero 
is  typically  a smii;  I ' i_aiit  i ,'iilai  m Cfi.j.  tin,_  tiicsc  values. 

Tile  lij'ic  lit  i/ciju ; Cl  ..  t‘t'  .1  i f Jt-i’t-'ii  l iavuis  ol 

luiictii'n  lUiii  gi'-uJicut  i 11 1 uriua  t i ai,  btcau.st;  t.lu  I'.'.i.i  i allowed  to  seJect 
eir.her  a s tep- 1 erigLh  algoriiiim  based  on  u.sing  only  iiiiiclicii  vaiue.s,  or  a 
step- J (.  ngtli  algorlLlaii  that  requires  iioth  l.i'ie  iunc  ti  uii  and  gradienl  (see 
Sei'Lion  The  c.alL  lo  St'i'N  must  tiierefore  allow  die  return  of  the 

follovv’ing  combi  nat  inns : both  tne  iimeuion  and  giadien:,  only  the  func- 

tion valun^  or  oniy  the  giadiem,  I'he  required  flexibility  is  supplied 
by  an  integer  paraiHeter  i IFIAG)  in  the  calling  sequence  of  3FUN,  which 
is  set  by  QNMDER  b.ifore  each  cllII  to  SPUN.  in  particular,  IFLAG-TJ  means 
that  only  the  function  value  is  to  be  computed,  IFLAG-1  moans  that  only 
the  gradient  vector  is  computed,  and  T.FIA&-S  means  that  both  the  func- 
tion and  gradient  should  be  returned.  Tlie  user  of  QNMDER  must  accordingly 
provide  a subroutine  to  serve,  as  the  actual  jiarameuer  corresponding  to 
SPUN,  such  that  the  correct  logic  is  followed  for  each  value  of  IFLAG. 

Hov/evor,  it  was  felt  that  such  flexibility  is  inappropriate  for  the 
user  of  the  ''easy-to-use''  routines,  who  should  not  be  required  to  program 
anything  but  an  absolutely  straightforward  routiae.  Consequently,  the 
’easy-to-use"  routine  that  calls  QElilDER  illustrates  a different  relation- 
ship with  the  user- supp  li  nd  subroutines.  Xlie  "e44sy-L1.-u.se"  subroutine 
always  chooses  the  .step  length  algorithm  that  requires  both  function  and 
gradient  values;  thus  IFlAG-i  in  all  calls  to  SFlUl,  so  tl;at  IFIAG  is 
redundant.  Requiring  tlu.  user  of  the  "easy-to-use"  routine  to  include 
a useless  j'aranieter  in  rhe  calling  sequence  of  nis  proldem  function  sub- 
routine iatrocluce.s  undesirable  complications,  on  the  ther  hand,  the 
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subroutine  QNMDER  should  not  be  altered  simply  to  make  it  compatible  with 
the  'easy-to-use''  version.  Hence,  a library  subroutine,  named  SFUN;1,  is 
provided  to  serve  as  the  parameter  corresponding  to  SI’UK  in  the  "easy- to 
use"  call  of  QNMDER.  l1ie  subroutine  SFUN2  then  calls  the  user- supplied 
routine,  which  now  must  have  not  only  a fixed  calling  sequence,  but  also 
a designated  name.  FUN...  The  ■'.ailing  sequence  of  FUNl;  is  then  of  the 
minimal  form  FHN2( N, X, F, G)  , where  N is  the  number  of  variables,  X is 
the  vector  of  variables,  F is  the  function  value  and  G the  gradient 
vector  evaluated  at  X . 

Similar  conventions  are  observed  throughout  the  library  in  order  to 
communicate  with  the  user.  The  existence  of  two  levels  at  which  the  user's 
routines  are  called  ~ directly  (with  fixed  calling  sequence)  and  indirectly 
(with  fixed  calling  sequence  and  designated  names)  — allows  both  flexi- 
bility in  the  general  case  and  simplicity  in  the  "easy-to-use"  case.  The 
requirement  of  a designated  name  .n  the  "simple"  case  makes  it  less  con- 
venient for  the  user  to  call  the  "easy-to-use"  library  routines  to  mini- 
mize different  functions  during  the  same  run,  but  this  sacrifice  was 
believed  to  be  insignificant  in  light  of  the  extreme  complication  of  pro- 
viding "easy-to-use"  routines  without  such  a restriction. 

3.9  Documentation 

The  importance  of  good  documentation  is  always  stressed  in  discussions 
of  mathematical  software  (see,  for  example,  Shampine,  Watts  and  Davenport, 
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iHiciiiiictr.  it  i loi  proSc-ut  libraty  rouLiac;.  occurs  in  two  torms. 

Kxtoindi  Jticu.iientai.  i o:i  prt)viclus  lIk  nsc-r  of  an  alv.orithn  with  in  fori.iat  i on 
r.c'v  i.-ssary  for  r.ho  1'fii.iniu  ni  an!  'lao  ol  t.l;^  j .id  i vi  cinal  rouLinu.s. 

In-  'inn  doc.unu  nlat  ioi:  apin-aiK  .•.■iiiiin  the  .-,ource  lent  uj  tiiC  prugraiii  tind  is 
intended  to  euabh*  p vogran  riU;r;'  wta.  a.ay  not  tiave  been  involved  in  Itie  original 
de.'.iua  CO  modify  tiv  rcutine  in  light  of  any  nev.  developmenta. 

As  a ! r.  ady  explained  in  Section  the  .-;t  run  are  of  optimization 

problems  reqeireb  that  the  iib..ary  contain  a variety  oi  routities.  Cite  user 
must  utilize  his  knovieijge  oi  tliv  problem  i e.  g.  , the  objective  lunction  is 
a snin  OI  squares),  aiid  the  ainount  and  type  of  at  ailable  i nformation  concern- 
ing the  problem  functions  fe.g.,  analytic  first  derivatives  can  be  computed) 
to  Select  the  aiproiii'iate  routine.  Unfortunately,  to  the  inexperienced  user 
it  n.ay  afipear  that  there  are  a bevildevirig  number  of  apparently  similar  rou- 
tines. in  order  to  assist  the  aser  in  selecting  the  aiJpropriate  routine, 
algorithm  selection  charts  are  provided  wiiii  the  libi:ary;  this  type  of 
decision  procedure  occurs  in  program  libraries  tor  solving  various  sorts 
f ptidilems  (for  e,--.ample,  the  EISP.AUK  collection),  the  selection  chart 
illustrated  in  Figure  ! concerns  tin.-  uncons  trained  oprlmization  of  a func- 
*ion  1-1  whoE'j  vector  of  first  derivatives  is  g(x)  and  whose  matrix  of 
secomi  derivatives  is  . 'ihe  user  mutt  decide  vihrtlicr  to  choose  a 

r.HUinc  tha’  computes  [■':•.),  _cu"  F :)  and  g(x)  , _or  l(x),  g(x)  and 
C x)  . ihe  basic  priticiple  is  to  utilize  as  much  derivative  information 
as  can  be  efficiintiy  computed,  because  iiietliods  using  the  higher  deriva- 
tives arc  more  robust  and  give  more  xiuonnation  about  the  quality  of  the 
solution  obtained.  Each  patli  th.rough  the  chart  leads  to  the  name  of  the 
.iihrouriues  tliat  dioul<l  be  used. 
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.ECTION  CHART  : UNCONSTRAINED  MINIM  I ZAT  ION 


Kacli  library  subroutine  is  described  in  detail  hy  ait  individual  report 
(sec,  fur  example,  OLll  et  al.,  I'.t'i'ba',  in  .diich  an  effort  lias  been  made  to 
provide  an  explanation  of  ilie  flicoretical  algoritluiis  as  well  as  the  com- 
puter implement  at  ions.  The  in-line  documeniation  includes  a block  of 
exp lanatiiry  comments  at  the  heginnirij/  of  every  subroutine,  containing  a 
discussion  of  the  numericai  procedure,  references^  and  a description  of 
eacii  p,iraiiie ter.  fa  addition,  comments  surround  every  non-obvious  block 
o f cudt . 

'ilic  need  to  extend  or  modify  library  routines  is  so  [irevelant  that 
it  must  be  included  as  a consideration  from  the  outset.  A significant 
aid  to  extension  (and  portability)  is  the  observance  of  rigid  textual  and 
programming  conventions  concerning  declarations,  the  ordering  of  formal 
parameters,  indentation,  speci fication  of  comments,  placement  of  machine- 
or  precision-dependent  quantities,  etc.  Such  policies  alsr  allow  a text 
editor  to  make  straightforward  changes  to  alter  the  code  for  various 
machines  and  precision  ranges.  Tne  rules  of  indentation  allow  a visual 
perception  of  program  structure,  and  assist  modification  by  indicating  the 
control  block  level  of  every  statement.  The  appearance  of  a library  pro- 
gram should  reflect  the  same  principles  of  clear  structure  that  guide  the 


analysis  of  riit;  methods  and  creation  of  the  implementations. 


A DotaiieJ  Illustration  — InipIcTUMitaLiun  of  a Step- l.euj^th  Algorithm 


li.  1.  Need  for  a s top- length  al^orltlim 

A step- length  algorithm  can  be  defined  in  general  terms  as  a pro- 
cedure that  solves  the  following  problem;  given  an  initial  point,  x , 
and  a direction,  p , find  a positive  scalar  step  u along  p , such 
that  at  X i >p  certain  conditions  are  satisfied  with  respect  to  those 
conditions  at  x . In  all  the  cases  to  be  considered  here,  the  desired 
step  can  be  viewed  as  an  approximation  in  some  sense  to  a , a local 
minimum  of  a specified  function,  say  0(x)  , along  the  given  direction. 

A step- length  algorithm  is  involved  as  a sub-step  in  almost  all 
unconstrained  optimization  algorithms.  This  point  is  emphasized  be- 
cause there  exists  a misunderstanding  that  a step- length  algorithm,  or 

''line  search",  is  equivalent  to  a procedure  to  locate  a highly  accurate 
* 

estimate  of  a for  the  given  ^ . For  example,  an  advantage  of 
Davidon's  "optimally  conditioned"  quasi-Newton  algorithm  Davidon, 

1975)  i-S  said  to  be  that  it  requires  "no  line  searches"  (Hillstrom,  19 1'' 
Nazareth,  1976).  However,  any  method  that  does  not  alter  the  search 
direction  during  an  iteration,  yet  requires  a decrease  in  some  ^^(x)  at 
each  iteration,  must  contain  a step- length  algorithm;  such  a method 
could  be  described  as  "without  line  searches"  only  if  it  required  otie 
function/gradient  e' aluation  per  iteration  without  exception. 

A step- length  algorithm  appears  not  only  in  methods  for  general 
unconstrained  optimization,  but  in  many  other  applications  as  well, 

Ih  linearly  constrained  problems,  the  step  length  is  frequently  chosen 
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either  to  satisfy  some  termination  criteria  with  respect  to  the  objective 
function,  or  as  the  maxinuim  feasible  step.  For  nonlinearly  constrained 
problems,  the  termination  of  the  step- length  algorithm  may  be  based  on 
some  modified  function  that  involves  the  original  objective  and  con- 
straint functions.  The  step- length  algorithms  in  the  library  should, 
therefore,  be  implemented  with  the  intention  that  they  will  be  required 
in  widely  varying  contexts. 

h.2  Choice  of  a step -length  algorithm 

In  this  section  we  note  the  results  in  the  library  design  of  a par- 
ticular decision  to  retain  flexibility. 

Good  step- length  algorithms  are  typically  based  on  safeguarded 
parabolic  or  cubic  approximations  to  the  behavior  of  the  given  function 
along  the  search  direction  (see  Gill  and  Murray,  197^a) ; the  cubic  pro- 
cedure utilizes  function  and  gradient  values  at  every  trial  point,  while 
the  parabolic  procedure  uses  only  function  values  except  possibly  at  the 
initial  point.  The  choice  of  a parabolic  or  cubic  method  is  sometimes 
straightforward.  For  example,  within  a Newton- type  algorithm  based  on 
exact  or  differenced  second  derivatives,  the  slightly  more  robust  cubic 
step-length  algorithm  is  always  used,  since  the  cost  of  evaluating  the 
gradient  at  every  trial  point  in  addition  to  the  function  is  assumed  to 
be  reasonable.  On  the  other  hand,  within  an  algorithm  based  on  finite- 
difference  approximations  to  first  derivatives,  the  parabolic  algorithm 
is  chosen  because  of  the  excessive  cost  of  the  n function  evaluations 
required  to  estimate  the  gradient  as  the  step-length  algorithm  proceeds. 
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However,  with  a quasi-Newton  method  tliat  uses  analytic  gradients, 
tha  best  choice  of  step-length  algorithm  (''best''  in  the  sense  of  'with 
fastest  execution  ; ime")  depends  on  the  expeme  involved  in  computation 
of  the  gradient  vector  relative  to  a calculation  of  the  function  value, 
since  the  parabolic  method  will  normally  lead  to  faster  execution  time  if 
c?  Icuiation  of  the  gradient  is  significantly  more  expensive  than  evalua- 
tion of  the  function.  When  developing  the  library  routines  for  this  case, 
it  was  felt  that  the  choice  of  step-length  algorithm  should  not  be  made 
a priori  for  the  knowledgeable  user  (in  contrast  to  the  ''easy-to-use"  case; 
see  Section  3*6). 

The  inclusion  of  this  flexibility  has  several  implications  with  regard 
to  the  programming  detail  within  the  relevant  routines.  First,  the  name  of 
a step-length  subroutine  must  be  included  as  a formal  parameter  in  the 
calling  sequence  of  the  library  subroutines  within  which  the  choice  is 
permitted.  Second,  because  either  method  can  be  chosen,  the  calling  se- 
quences of  the  subroutines  that  carry  out  both  the  parabolic  and  cubic 
procedures  must  be  the  same.  Finally,  special  logic  needs  to  be  included 
in  the  calling  subroutines  to  test  which  method  was  selected,  so  that  the 
necessary  gradients  can  be  computed  upon  termination  of  each  line  search 
if  the  parabolic  method  was  used.  Throughout  the  library,  there  are  simi- 
lar instances  where  non-obvious  consequences  follow  an  apparently  simple 
design  decision. 

1.3  Provision  of  user-defined  parameters 

The  calling  sequences  of  several  general  optimization  subroutines 
include  two  parameters  that  define  in  part  the  termination  criteria  of 
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the  i toi> -1  ength  algorithm  Hero  wi'  ceiitiick.'r  why  theac!  Lwu  specific 
levels  of  ILexibility  have  been  provided. 

Tlie  first  parameter,  | (E'IA),  0 ^ T|  • 1 , controls  the  accuracy 

to  which  the  step-length  procedure  locates  a local  minimum  of  the  given 
function  along  the  search  direction.  l.et  0 denote  the  function  to  be 
minimi.^ed,  and  let  g denote  its  projected  gradient  along  the  search 
direction.  For  the  cubic  case,  the  step-length  algorithm  normally 
terminates  with  i if: 


|g(x  rap)  i < ■,  |g(x)  ] , 


For  the  quadratic  case,  the  algorithm  normally  terminates  when  a local 
minimum  of  0 along  p is  known  to  be  contained  in  the  interval  [a,b] 
and 


0yx  +gp)  - 0(x  i ap) 


< 


|g(x) 


i.e.,  when  the  lincarixcd  appro;:iniation  to  the  projected  gradient  at 
X +ap  satisfies  the  same  test  as  g(x+ap)  in  tlie  cubic  case. 

The  computation  time  required  to  solve  a given  problem  will  vary 
with  the  choice  of  ■,  ; reducing  r,  tends  to  decrease  the  number  of 
iterations,  but  increase  the  number  of  evaluations  of  the  problem  func- 
tions, with  the  reverse  effect  from  increasing  i|  . The  "best"  value 
of  T]  (that  allows  the  problem  to  be  solved  in  the  shortest  time)  de- 
pends on  the  problem.  The  variation  in  the  optimal  is  small  for 
most  problems,  and  there  is  usually  little  decrease  in  efficiency  if 
-i  is  set  to  an  averaged  optimal  value  based  on  extensive  computational 


experimentation  (as  in  the  easy-to-use  routines).  However,  on  some  cate- 
gories of  problems,  the  optimal  value  of  may  differ  significantly 
from  this  average.  In  such  instances,  a considerable  improvement  results 
by  allowing  the  user  the  freedom  to  adjust  ij  to  suit  the  problem,  and 
consequently  the  parameter  ETA  is  retained  in  the  calling  sequence. 

The  second  parameter  to  be  considered,  which  strongly  influences 
the  performance  of  the  step-length  algorithm,  is  X (STEPMX)  , an  upper 
bound  on  the  norm  of  the  step  to  be  taken  during  any  iteration  (i.e.,  the 
final  a must  satisfy  i|o:p!|  < X).  There  are  several  reasons  for  includ- 
ing this  parameter  in  the  step-length  algorithm: 

( 1)  to  prevent  overflow  in  the  user-supplied  problem  function 
subroutine ; 

(2)  to  increase  efficiency  by  evaluating  the  problem  functions 
only  at  ''sensible"  values  of  x ; 

(5)  to  prevent  the  step-length  algorithm  from  returning  an 

inordinately  large  step  because  no  smaller  step  satisfies  the 
relevant  convergence  criteria.  It  is  not  always  appreciated 
that,  even  if  a function  is  unimodal  in  E*^  , there  may  exist 
directions  along  which  it  is  monotonical ly  decreasing; 

(X)  to  attempt  to  force  convergence  to  the  local  minimum  nearest 
to  the  initial  estimate. 

Because  the  user  should  be  in  the  best  position  to  specify  a good 
choice  of  X , this  parameter  is  also  included  in  the  outermost  calling 
sequence.  For  the  easy-to-use  routines,  however,  X is  automatically 
set  at  max{  10,  “x^^ ^ i| ) , where  x^^^ 


XX 


is  the  initial  value  of  x . This 


choice  o£  X is  dependent  on  the  prohii.-in  scaling,  and  hence  the  per- 
formance of  the  easy-to-use  routines  may  he  degraded  on  badly  scaled 
problems. 

i;.l4  Control  structure  for  a step-length  algorithm 

A step-length  algorithm  should  be  implemented  with  the  knowledge 
that  it  will  appear  almost  universally  throughout  the  library  (see 
Section  4.1).  Conceptually,  the  execution  of  a given  step-length  pro- 
cedure is  identical  in  any  context,  and  one  might  hope  to  implement  the 
parabolic  and  cubic  step- length  algorithms  as  invariant  subroutines", 
this  section  considers  the  complications  of  such  an  implementation, 
which  initially  seems  straightforward. 

One  significant  element  of  a step- length  algorithm  involves  the 
retention  and  communication  of  information.  In  the  unconstrained  case, 
the  values  of  the  objective  function  and  its  gradient  vector  at  the  best 
point  must  be  stored  as  the  linear  search  proceeds,  to  be  communicated  to 
the  outer  iteration  upon  termination.  Because  the  best  point  found  is 
not  necessarily  the  most  recent  point  at  which  the  function  was  evaluated, 
the  step-length  routine  also  requires  storage  to  contain  the  most  recent 
function  and  gradient  values.  In  this  context,  the  information  to  be 
retained  during  execution  of  each  linear  search  involves  only  the  func- 
tion upon  which  the  termination  criteria  of  the  step-length  algorithm  are 
based  (hereinafter  referred  to  as  the  "intermediate  function").  However, 
in  other  instances  the  information  associated  with  a trial  point  is  not 
restricted  to  the  value  and  gradient  of  the  intermediate  function. 
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example,  in  the  nonl 
retained;  in  nonlinearly  constrained  problems,  the  data  to  be  stored  include 
at  the  least  the  values  and  gradients  ol  the  original  objective  and  con- 
straint functions.  It  is  difficult  to  imagine  an  invariant  step-length 
subroutine  tliat  would  be  suitable  for  all  situations,  without  resort  to 
prt)gramming  contortions  involving  complex  storage  arrangements.  The  neces- 
sity of  communicating  varying  sets  of  information  eliminates  the  desira- 
bility of  developing  a universal  step-length  module. 

However,  it  is  possible  to  isolate  a second  element  of  a step-length 
algorithm  that  is  not  dependent  on  the  information  connected  with  each 
point;  this  aspect  is  the  calculation  of  the  next  trial  point,  which  is 
based  on  a safeguarded  polynomial  approximation  procedure,  and  utilizes 
quantities  relevant  only  to  the  intermediate  function.  Consequently,  the 
library  contains  fixed  subroutines  to  i.-nplement  this  second,  problem- 
independent,  portion  of  the  step-length  algorithms. 

The  fixed  subroutines  are  NEWPTC  (Gill  et  al.,  1976b)  and  NEWPTQ 
((lill  et  al.,  197'^  a-),  which  carry  out  a single  iteration  of  a step-length 
procedure;  NEWPTC  computes  the  next  trial  point  based  on  safeguarded  cubic 
approximation,  while  NEWPTQ  uses  safeguarded  quadratic  approximation.  The 
step-length  procedure  for  a given  application  involves  repeated  calls  to 
one  or  the  other  of  these  modules  by  an  outer  controlling  subroutine,  which 
is  tailored  to  the  particular  context;  the  outer  routine  retains  the  neces- 
sary information  as  the  iteration  proceeds,  and  provides  NEWPTC  and  NEWPTQ 
with  the  value  (and  possibly  projected  gradient)  of  the  intermediate  func- 
tion at  each  new  point.  In  this  way,  all  calls  to  user-supplied  or  library 


rovitiines  that  define  the  problem  functions  occur  outside  NEWPTC  and  NliMP'rQ. 
Such  a control  structure  was  suggested  by  the  "reverse  communication" 
mechanism  of  Lawson  and  Kvogh,  dis['layed  in  the  differential  equation 
solver  DVLX5  (see  Krogh,  ly6y). 

The  displeasing  feature  of  the  chosen  arrangement  is  that  a parameter 
of  NEWPTC  and  NEWPlf)  must  indicate  whether  the  current  call  is  the  first 
within  an  iteration,  and  whether  certain  other  special  conditions  hold. 

The  concept  of  control  defined  by  a co- routine  (see  Dahl  and  Hoare,  197^) 
would  be  a more  accurate  reflection  of  the  nature  of  the  step-length 
algorithm,  but  does  not  exist  in  current  Fortran. 

The  overwhelming  advantage  of  the  given  design  is  the  separation  of 
the  two  distinct  tasks  of  obtaining  the  next  trial  point  and  maintaining 
the  associated  information.  The  implementation  of  the  inner  logic  of  the 
step-length  algorithms  — safeguarded  cubic  or  parabolic  approximation— 
thus  is  independent  of  the  dimensionality  and  data  structures  of  the  par- 
ticular optimization  problem  within  which  a step-length  procedure  is 
required. 
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llie  design  principles  underlying  an  optimization  program  library  ' 

have  been  presented,  in  order  to  indicate  the  kinds  of  considerations 

and  compromises  involved  in  such  a major  software  project.  It  is  i 

impossible  to  achieve  a ''perfect''  library,  but  at  least  every  decision  i 

in  the  present  system  has  been  made  based  on  an  analysis  of  the  same  j 

set  of  aims.  I’here  will  inevitably  be  unforeseeable  future  developments 

and  modifications;  in  fact,  many  changes  have  already  occurred  since  the 

first  routines  were  produced  in  197^’  However,  it  is  hoped  that  the  care 

put  into  the  library  design  will  continue  to  allow  consistent  and  well- 

structured  extensions.  ' 

i 

i 

i 
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I . Appendix  ( an  example  of  documentation) 

SUBROUTINK  UCNDQl 

An  easy-to-use  quasi- Newton  algorithm  to  find 
the  unconstrained  minimum  of  a function  of 
N variables  using  function  values  only. 

Philip  E.  Gill,  Walter  Murray,  Susan  M.  Picken, 

Margaret  H.  Wright  and  Hazel  M.  Barber 

Ref.  No.  EVlO/O/Fortran/ll/Yp 

NPL  ALGORITHMS  LIBRARY  Fortran  SUBROUTINE  UCNDQl 

1.  Purpose 

To  minimize  a function  F( xl, x2, . . . , xN)  of  the  N independent 
variables  xl,x2, . . . ,xN,  using  a quasi- Newton  method.  The  sub- 
routine UCNDQl  is  intended  for  functions  which  are  continuous  and 
which  have  continuous  first  and  second  derivatives  (although  itvv’ill 
usually  work  if  the  derivatives  have  occasional  discontinuities). 

The  user  must  supply  an  initial  estimate  of  the  position  of  a 
m.inimum  and  a subroutine  which  will  calculate  F(x)  at  any  point 
X - ( .X 1 , x2, . . . , xN ) . It  is  essential  that  the  subroutine  supplied  by 
the  user  has  the  name  FUNl  and  has  the  correct  specification. 


SUBROUTINE  UCNIX^l  is  an  casy-Lo-use  version  of  SUBROUTINE  QNMDIF, 
NPL  AI.GORITHMS  1.1BR.VRY  Ref.  No.  El/Op/F  and  is  intended  for  users 
wlio  have  absolutely  no  knowledge  of  cither  optimization  or  the 
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behavior  of  their  problem.  It  mu.sL  be  emphasised  that  a minimization  i 

can  always  be  perfonned  more  efficiently  if  the  SUBROUTINE  QNMDII’  is 
used  and  the  parameters  arc  chosen  to  suit  the  problem  being  solved. 

llie  user  should  also  be  aware  of  the  alternatives  to  SUBROUTINE  i 

UCNOQl.  Tlie  SUBROUTINE  VCFVQd,  NPL  ALGORIIUMS  LIBRARY  Ref.  No. 

Eh/Oj/F,  is  available  which  requires  the  user  to  calculate  function 
and  gradient  values.  There  is  also  SUBROUTINE  UCFDN2,  NPL  ALGORITHMS 
LIBRARY  Ref.  No.  E^i/oS/F,  which  requires  the  user  to  calculate  func- 
tion and  gradient  values,  and  SUBROUTINE  UCSDN2,  NPL  ALGORITHMS  j 

LIBRARY  Ref.  No.  eH/OT/F,  a similar  routine,  except  that  it  also 
requires  the  user  to  calculate  second  derivatives. 

i' 

2.  Description 

From  a starting  point  supplied  by  the  user,  a sequence  of  points  • 

is  generated  which  is  intended  to  converge  to  a local  minimum  of 

F(x).  These  points  are  generated  using  estimates  of  the  gradient 

and  curvature  of  the  objective  function.  Finally,  an  attempt  is 

made  to  verify  that  the  final  point  is  a minimum. 

Specification 

SUBROUTINE  UCNDQ1(N,  X,  F,  IFAIL,  W,  LW) 

INTEGER  N,  IFAIL,  LW 

REAL  F I 

REAL  Xf  N) , W( LW)  ] 

! 

i 


^0 


Parametors 


Input  parameter 

N - integer,  a cwnstant  w/ii  h ny.iat  in-  sei  tm  entry  to  the 

number  oI  independent  variahles  in  the  funetjon  to  be 
minimized. 

Input  and  euiiput  parameter 

X - RI'AI,  array,  of  length  N. 

Oil  entry  to  lICNIX^l,  X(j)  muxc  liave  been  ret  by  tlie  user  to 
an  estimate  of  tiie  JLh  corapunenl.  of  the  position  oi  the 
minimum  ij  - 

After  a normal  exit,  X(j)  contains  the  Jth  component  of 
the  position  of  the  mininium. 

Output  parameters 

F - REG\h.  On  exit,  F gives  the  value  of  F'x)  corresponding 

to  the  final  point  stored  in  X. 

IFAII.  - INTEGER.  On  exit,  TFAIL  = 0 indicates  a successful  call 
of  the  subroutine  UCNIXIl.  For  other  exit  values  of  IFAIL 
and  tlieir  meanings  see  Error  Indicators.  N.  B.  Tlie  user 
must  test  the  value  of  IFAIL  to  determine  whether  the 
subroutine  lias  run  successfully . 

Workspace  parameters 

W - REAl.  array  of  length  L.W,  used  as  workspace  by  UCNIXJl. 

Elements  of  W used  by  UCNIXJI  are  W(j),  .1  1,  1 •N>N‘'  N- 1 ) 

or  J%] ,2, . . . , 1 1 if  N_1 . 


.1 
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- TNTKli'KR,  a constant  whicli  gives  the  actual  lengtli  of  the 
array  W as  declared  in  the  program  unit  from  which  UCNUQl 


is  called.  This  number  must  be  at  least 
or  n.  if  N 1. 


iO'N  -I-  N'f  N-1)/.. 


'j.  User  supplied  subroutine 

I’l'Nl  - It  is  essential  that  the  subroutine  supplied  by  the  user 
to  calculate  the  value  of  the  function  has  the  name  FUNl 
and  has  the  specification:- 
SUBROUTINE  FUN1(N,  X,  F) 

INTEGER  N 
REAh  F 
REAL  X(N) 

The  function  value  at  the  point  defined  by  the  array  X(l), 

1 = 1^2^..., N must  be  stored  in  F. 

The  elements  of  the  array  X and  the  integer  N must  not  be 
altered  in  ITJNl. 

' • Error  Indicators 

IFAIl.  1,  parameter  outside  expected  range.  This  failure  will 
occur  if,  on  entry,  NCl,  LW-"(10‘N  f N'(N-1)/2)  or  LW-^11  if 
N 1. 

IFAIL  there  have  been  ltOO''<N  function  evaluations,  yet  the 

algoritlim  does  not  seem  to  be  converging.  The  calculations 
can  be  restarted  from  the  final  point  held  in  X.  The  error 
may  also  indicate  tliat  F(x)  has  no  minimum. 


\ 
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IFAIL  - 3»  the.  condiLiouii  fur  a iiinimum  have  not  ail  been  met 

but  a lower  point  could  nut  be  found  and  the  rllini!T^i^ation 
has  failed. 

IFAfL  = kf  an  overflow  has  occurred  during  the  computation.  Tiiis 
is  an  unlikely  failure,  but  if  it  occurs  tiie  user  sliould 
restart  at  the  latest  point  given  in  X. 

IFAIL  = p,  the  conditions  for  a minimum  liave  not  all  been  met 
but  the  minimization  has  probably  worked. 

IFAIL  - o,  the  conditions  for  a minimum  have  not  all  been  met 
but  tlie  minimization  has  possibly  worked. 

IFAIL  = , the  conditions  for  a minimum  have  not  all  been  met 

and  the  minimization  is  unlikely  to  have  worked, 

IFAIL  - 8,  the  conditions  for  a minimum  have  not  all  been  met 
and  the  minimization  is  very  unlikely  to  have  worked. 

If  the  user  is  unsatisfied  with  the  solution  it  is  worth  restarting 
the  calculations  from  a different  starting  point  (not  the  point  at 
which  the  failure  occurred)  in  order  to  avoid  the  region  which 
caused  the  failure.  If  persistent  trouble  occurs  and  the  gradient 
can  be  calculated  it  may  be  advisable  to  use  subroutine  UCFIX)2  or 
UCFDNi:  (see  Purpose). 

7.  Auxiliary  Routines 

Routines  named  SFb'Nl,  SPRNTl,  DOT,  LDLTSL,  MDCOND,  NMDCHL,  MODQIL, 

QNMDIF,  iVPPGRD,  FRMDEL  and  [JNSCll  are  called  by  UCNDQl.  Auxiliary 
routine  NEWPTQ  is  called  by  LINSCH. 

The  texts  of  the  first  two,  together  with  UCNIX)1,  are  given  in 


this  document.  The  texts  of  the  others  are  given  in  NPL  ALGORITHMS 
LIBRARY  documents,  reference  numbers  eJ|/02/F,  E4/O5/F,  EL/06/F  and 
E-V/15/F. 


8.  storage 

There  are  no  internally  declared  arrays. 

9.  Timing 

The  number  of  iterations  required  depends  on  the  number  of  variables, 
the  behaviour  of  F(x)  and  the  distance  of  the  starting  point  from 
the  solution.  The  number  of  operations  performed  in  an  iteration  of 
UCNDQl  is  roughly  proportional  to  N**2.  In  addition,  each  iteration 
makes  at  least  N+1  calls  of  FUNl.  So,  unless  F(x)  can  be  evaluated 
very  quickly,  the  run  time  will  be  dominated  by  the  time  spent  in 
FUNl. 

10 . Accuracy 

When  a successful  exit  is  made  then,  for  a computer  with  a 
wordlength  of  t decimals,  one  would  expect  to  get  about  t/2  - 1 
decimals  accuracy  in  x and  about  t - 1 decimals  accuracy  in  F, 
provided  the  problem  is  reasonably  well  scaled. 

11.  Further  Comments 

Ideally,  the  problem  should  be  scaled  so  that  tiie  minimum  value 
of  F(x)  is  in  the  range  (-1,  +1),  and  so  that  at  points  a unit 
distance  away  from  the  solution  F is  approximately  a unit  value 
greater  than  at  the  minimum.  It  is  unlikely  that  the  user  will  be 
able  to  follow  these  recommendations  very  closely,  but  it  is  worth 
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trying  (by  guesswork),  as  sensible  scaling  will  reduce  the  difficulty 
of  the  minimization  problem,  so  that  UCNDQl  will  take  less  computer 
L ime . 

12.  Ke'vnjords 

Easy-to-use,  function- only  method,  mathematical  programming,  maximi- 
zation, minimization,  optimization,  quasi-Newton  method,  SUBROUTINE 
QNMDIF,  unconstrained. 

IMPLEMENTATION  INFORMATION 

Tliis  program  (except  for  the  first  line)  has  been  written  in  ANSI  Fortran 
and  has  been  tested  on  a CDC6>00  computer. 

The  machine- dependent  constant  EPSMCH  must  be  set  in  the  first  instruction 
in  SUBROUTINE  UCNDQl,  in  SUBROUTINE  QNMDIF,  in  SUBROUTINE  APPGRD,  in 
SUBROUTINE  NMDCHL,  and  in  SUBROUTINE  FRMDEL.  The  machine- dependent  constant 
RMAX  must  be  set  at  the  beginning  of  SUBROUTINE  QNMDIF  and  of  SUBROUTINE 
MODCHL.  The  constant  EPSMCH  is  the  smallest  positive  real  number  such 
that  1.0  + EPSMCH  > 1.0,  and  RMAX  is  the  largest  positive  real  number 
such  that  both  RMAX  and  -RMAX  can  be  held  in  the  machine.  The  example 
program  has  set  EPSMCH  = 2.0Ef0*'’(-4'7),  and  RMAX  = l.OE+294  which  are 
ttie  correct  values  for  the  CDC  65OO  computer. 

EXAMPLE 

Minimize  the  function 

F(xl,x2,x3,xi4)  = (xl  -i  10.0^x2)*<*2  + ‘;.0^‘(x3  - x4)**2  + ( x2  - 2.0*x3)-'‘k 
+ 10.0*(xl  - xk)**k, 

starting  from  the  initial  estimate  X = (3.O,  -1.0,  0.0,  l.O). 


This  problem  has  the  solution  F(x)  = 0.0  at  x = (O.O,  0.0,  0.0,  O.O), 
and  is  rather  difficult  in  that  the  problem  has  a long  valley  along 
which  the  function  is  slowly  varying. 


PROCEDURES  AND  EXAMPLE  PROGRAM  TEXT 

Now  follows  the  text  of 

a)  the  program  for  the  example  problem, 

b)  a user  specified  subroutine  FUNl,  calculating  the  function  value  of 
the  example  problem, 

c)  the  subroutine  UCNDQl, 

d)  two  auxiliary  subroutines  SFUNl  and  SPRNTl. 

The  other  auxiliary  routines  required  may  be  found  in  NPL  ALGORITHMS 
LIBRARY  documents  reference  numbers  Eh/02jF,  Eh/OJ/jF,  E1|/06/f  and 
Eh/15 /F. 

This  program  has  been  run  successfully  on  a CDC  65OO  computer. 
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